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SUMMARY 


A finite element hydroelastic analysis formulation is developed on the 
basis of Toupin's complementary variational principle. Emphasis is placed on 
the special case of an incompressible fluid model which is applicable to pro- 
pellant tank hydroelastic analysis. A concise fluid inertia representation re- 
sults from the assumption of incompressibility and the hydroelastic equations 
reduce to a simplified form associated with the structure alone. The efficien- 
cy of the incompressible hydroelastic formulation is enhanced for both fluid 
and structure by introduction of harmonic reduction as an alternative to Guyan 
reduction. The theoretical developments are impelemented in the NASTRAN 
Program and the technique is verified and demonstrated as an efficient and 
accurate approach with a series of illustrative problems including the 1/8- 
scale Space Shuttle external tank. 


INTRODUCTION 


The increasing complexity of launch vehicle configurations, particularly 
in the case of the Space Shuttle, recently has stimulated considerable interest 
in the dynamic behavior of liquid-filled tanks. The prevention of coupled 
structure -propulsion instability (pogo), for example, requires very complete 
and accurate mathematical models for the calculation of propellant tank hy- 
droelastic modes (ref. 1) in the frequency range of concern (2-50 Hz for the 
Space Shuttle ). 

Various fluid modeling techniques for hydroelastic analysis (refs. 2-7) 
have been mechanized for digital computation. These techniques range from 
finite element and finite difference techniques to approximate analytical 
approaches taking advantage of the properties of the fluid velocity potential 
and the consequences of Green's theorem. Most hydroelastic analysis methods , 
however, are limited to special geometric configurations. Furthermore, al- 
though theoretically rigorous, most contain deficiencies in computational 
economy and/or numerical accuracy. 
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In particular , the hydroelastic analysis technique used in NASTRAN 

7) is deficient in computational economy* This technique employs a 
large unsymmetrical eigenvalue problem formulated in terms of mixed fluid 
pressure and structural displacement generalized coordinates* In this for- 
mulation, the fluid coefficient matrices, derived by a Galerkin-type approach, 
are interpreted according to a structural analogy* The resulting fluid pseudo - 
mass and pseudo -stiffnes s matrices are recognized herein as flexibility and 
inverse mass (inertance) matrices, respectively, on the basis of Toupin*s 
complementary variational principle (ref. 8), This revised interpretation is 
central to the formulation of the hydroelastic problem presented here. Al- 
though the effort in the present work was initially directed towards alleviation 
of difficulties encountered in the NASTRAN hydroelastic analysis, it has been 
found that Toupin^s principle provides fundamental physical insight for a vari- 
ational approach to hydroelasticity and it provides a rigorous basis for the 
development of fluid finite elements * 

A derivation of fluid matrix equations on the basis of Toupin*s principle, 
and manipulation of the interacting structure equations to a form consistent 
with the complementary principle, \altimately results in a symmetric kinetic 
formulation for compressible hydroelasticity. A detailed development of the 
compressible formulation is not presented in this paper. The special case of 
incompressible hydroelasticity is of primary interest as it is particularly 
applicable in the study of propellant tank dynamics. 

In the special case of incompressibility, the interior fluid pressure 
fluctuations are algebraically dependent on surface pressure fluctuations. Al- 
though this property is illustrated from the viewpoint of a matrix approach, it 
is well known as a consequence of the application of Green's theorem to a 
continuum potential description. Moreover, the surface pressure fluctuations 
are algebraically related to the fluid accelerations normal to the bounding sur- 
face, and this is equivalent to a statement of the '’surface flux" boundary con- 
dition. In addition, it is recognized on physical as well as mathematical 
grounds that there exists an over specification in the surface pres sure /accel- 
eration relationships which implicitly defines a compatibility condition* In 
very simple terms the condition which requires that under uniform pressure 
the fluid will not move manifests itself as a necessary singularity in the fluid 
inertance matrix. Upon elimination of the surface pressure singularity by 
introduction of pressure deviation dynamic variables, a concise fluid mass 
matrix is formed in terms of bounding- surface displacements alone; this 
matrix is simply added to the structural mass matrix, resulting in a symmet- 
ric kinematic set of hydroelastic dynamic equations. This description repre- 
sents a drastic reduction in system variables. 

It is noted in this paper that the sequencing of kinematic reduction oper- 
ations (e.g,, Guyan, etc* ) on fluid and structure is crucial to the economy of 
the analysis approach. In addition, a matrix harmonic reduction scheme is 
introduced as an efficient alternative to Guyan reduction for geometrically, 
but not nece ssarily structurally, axisymmetric configurations to further re- 
duce the number of system variables, 

A series of illustrative problems is presented to demonstrate the pre- 
sent incompressible hydroelastic analysis, which has been implemented in 
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NASTRAN. In particular the problems illustrate the use of harmonic reduc- 
tion, and the accuracy and efficiency of the hydroelastic formulation in gen- 
eral, by comparison with exact analytical and available test results. More- 
over, computation times for the standard NASTRAN and the present hydro - 
elastic analysis techniques are compared. 
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SYMBOLS 

generalized area matrix 

I 

fluid bulk modiolus 
flexibility matrix 
fluid, structural flexibility matrix 
structural elastic modulus 
internal structural generalized force 
matrix defined in equation (29b) 
identity matrix 
stiffness matrix 

structural stiffness matrix, constrained structural stiffness 
matrix [eq, (32d)] 

inertance matrix 

fluid, structural inertance matrix 
fluid mass matrix 

structural mass matrix, constrained structural mass matrix 
[eq. (32c)] 

cylindrical-shell bending -moment resultants [eq. (40a)] 
cylindrical -shell membrane -stress resultants [eq. (40a)] 
pressure, pressure deviation array 
hemisphere or cylindrical shell radius 
surface area 

complementary kinetic, potential energy function 
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reduction transformation [eq. (34)], harmonic transformation 
[eq. (36)] 

displacement array- 
volume 

complementary virtual work function 
shell thickness 

cylindrical shell axial dimension 
meridional, circumferential wave index 
surface outward normal unit vector 
pressure 

static pressurization level (gage) 

harmonic distribution pressure amplitudes 

[eq. (37)] 

radial coordinate in cylindrical reference frame 
time 

displacement, displacement vector 
surface displacement 

axial coordinate in cylindrical reference frame 
matrix defined in equation (28c) 
nondimens ional frequency 

stiffness constant for hemisphere [eq. (38)] 

Poisson's ratio 

fluid, structural density 

circular frequency 
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Operators 


d( ) 

total differential 

V* ( ) 

divergence 

V( ) 

gradient 

M ) 

partial derivative 

6( ) 

variation 

( ) 

total impulse / ( ) dt 


-CO 

(•) 

time derivative, d( )/dt 

Subscripts 


i,j»k 

denote elements of a matrix or vector 

Abbreviations 


DOF 

degrees of freedom 

CPU 

central processing unit 


THEORETICAL DEVELOPMENT 


The class of problems considered here consists of the interaction of 
irrotational , inviscid, compressible fluids with flexible structures for which 
both fluid and structural motions are assumed small compared to overall 
dimensions. The approach used to describe the dynamics of the fluid is a ^ 
finite element technique which utilizes energy expressions based on Toupin's 
principle (ref. 8). A detailed derivation of Toupin's principle and a comple- 
mentary form of Hamilton's principle are presented in reference 9 as conse- 
quences of a postulated complementary D'Alembert principle. 

The equation of motion of a fluid particle is 


and the constitutive relationship for an inviscid, compressible fluid is 


P = 


- BV • u 


( 2 ) 


181 


where v • u represents the dilatational strain. In order to obtain a fluid 
velocity expression, equation (1) is integrated to yield 


u = 



where p is the pressure impulse 


A 


P 



or 



( 3 ) 


(4) 


Complementary kinetic and strain energy (T and U , respectively) may now 
be expressed in terms of impulsive pressur^ (the ccfmplementary dynamic 
variable) as 


• VP) dV 

V 


(5a) 


(5b) 


The motion dependent and impulse dependent energy expressions are generally 
not equivalent; they are equivalent, however, for linear systems (ref. 10). 

The complementary virtual work performed by boundary surface displace- 
ments,“u'% is 


6W 

c 



6p (u^'* n) dS 


( 6 ) 


The concept of complementary virtual work (ref. 10) may be viewed as a con- 
sequence of a complementary D'Alembert principle (ref. 9). 

The above energy expressions substituted into the complementary form 
of Hamilton's principle 
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U,)dt 




6W dt 
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(7) 
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ultimately yields, after use of Green* s theorem and integration by parts, 
integrated expressions corresponding to the wave equation and natural boun- 
dary conditions* 

The complementary formulation is presently applied, however, as an 
approximate analysis tool. Consider an approximation of a fluid pressure 
(impulse) state expressed as a linear function of a finite set of generalized 
impulses* Let us also require that any such approximation contains spatially 
uniform pressure as a state describable by the chosen generalized impulses; 
this is analogous to the requirement in kinematic finite element analysis that 
assumed displacement states must contain rigid-body motions. The fluid 
complementary kinetic and strain energies resulting from the assumed pres- 
sure impulse states are therefore the quadratic functions 
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(8a) 


(8b) 


with the elements of the symmetric inertance matrix L and flexibility matrix 
C defined as 


L.. 






C.. = 


c 


(9a, b) 


The Ly are proportional to 1/Pf and the C^j are proportional to 1/B. The 
complementary virtual work is expressed as 
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dS 6p. 


( 10 ) 


For the special case in which the surface displacements are physically dis 
cretized, the complementary virtual work may be expressed as 


6W 
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k i 


(11a) 
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with an element of the generalized area matrix A defined as 


A 


ik 




dS 


(11b) 


Substitution of equations (8) and (11) into equation (7), with the appro- 
priate integrations by parts, results in the complementary Euler -Lagrange 
equations 



*ik“k 


( 12 ) 


By taking the time derivative of this expression noting equation (4), the Euler - 
Lagrange equations become 

j k 


This is the form of the fluid dynamic finite element equations for indi- 
vidual elements and stacked systems of elements. In the case of a stacked 
system of elements, the matrix A represents only bounding surface general- 
ized areas and u^ represents discrete surface displacements. The pressures 
Pj comprise the set of boundary surface and internal pressures; therefore the 
matrix A is rectangular. 

The above set of fluid equations is derived for NASTRAN (ref. 7) by a 
Galer kin-type approach which "constructs" a minimal principle. In this lat- 
ter approach, the matrix expressions are interpreted according to a mathe- 
matical analogy and physical insight is lost; that is, L^j is mathematically 
analogous to kinematic stiffness, Cjj is analogous to kinematic mass and the 
right-hand side is analogous to a kiri^ematic generalized force vector. 


A Symmetric Kinetic Formulation for 
Compressible Hydroelasticity 


In the general case of compressible fluid/structure interaction, the fluid 
is described in terms of the complementary form 


L^P + C^P = - a'^U 


(14a) 
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and the structure under fluid pressure excitation is described in the standard 
kinematic form 


M U + K U = AP 
s s 


(14b) 


Direct coupling of the above set in terms of the mixed pressure and structural 
displacement variables [as in the standard NASTRAN formulation (eq. (7))] 
results in an xinsymmetric eigenvalue problem. 

A symmetric formulation can be derived by the complementary principle 
or by an equivalent manipulation of the structural dynamic equations to the 
complementary form. If the latter approach is taken, internal structural 
generalized forces, Fg, are related to the structural displacements, U, ac- 
cording to 


K U = F (15) 

s s 

Suppose that Kg represents a stiffness matrix of an interacting supported struc- 
ture such that rigid body deflections do not occur; the transformation to inter- 
nal forces is defined as 


U . (16) 

Incorporation of this transformation into the hydroelastic equation set, equa- 
tion (14), yields the symmetric set of hydroelastic equations in terms of force - 
type variables. 
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(17) 


with the structural inertance and compliance matrices defined, respectively, as 

L = M‘^, C = K'^ (18) 

s s ' s s ' ' 


The formulation presented here provides a symmetric kinetic formula- 
tion for inviscid, compressible fluid hydroelastic problems for which efficient 
eigenvalue analysis techniques are applicable. Further discussion of the com- 
pressible formulation is not presently warranted since the class of problems 
of interest is limited to incompressible fluids interacting with flexible struc- 
tures. The alternate simplified kinematic formulation to be derived below is 
appropriate for this case. 
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A Symmetric Kinematic Formulation for 
Incompressible Hydroelasticity 


In the special case for which fluid compressibility is negligible (C£ ~ 0) 
the complementary fluid dynamic equations, equation (14a), reduce to a set of 
algebraic equations relating fluid pressures and boundary -surface accelera- 
tions. One is strongly motivated to utilize this quality to solve for the pres- 
sures in terms of surface accelerations --and ultimately obtain a kinematic 
formulation which is identical in form to a standard set of structural dynamic 
equations. 

Two fundamental properties of the fluid are recognized as consequences 
of the simplifying assumption of fluid incompressibility. These properties are 
necessarily inherent in both the continuum and the present matrix finite ele- 
ment descriptions. The first property requires that for an incompressible 
fluid the interior pressures are related to the surface pressures in a purely 
geometric sense; from the continuum viewpoint, this is a consequence of 
Green’s theorem applied to an incompressible fluid. The second property re- 
quires that the net flux (or normal flow) out of the fluid volume is zero. This 
latter property amounts to a statement of constraint on flow normal to the 
bounding surface, and simultaneously it is recognized as a compatibility con- 
straint on surface pressures. A simple physical statement of this second 
property is that, under the uniform surface pressure (impulse) state, the fluid 
surface will not deform and consequently, the fluid volume will not move. This 
point is illustrated mathematically by noting that in the uniform pressure state 
the pressure gradient, vp, is null throughout the fluid volume; hence the com- 
plementary kinetic energy, equation (5a), is null. 

The complementary fluid matrix equation set for the special case of an 
incompressible fluid [see equation (14a) for Cf = 0] in a conveniently partitioned 
form is expressed as 



The pressure partitions, P£, and P^ correspond to a single surface 
reference pressure or free-surface pressure set, the complement of the total 
surface pressure set, and the internal fluid pressure set, respectively; the 
displacement partitions IJ£ and Ug correspond to a single surface reference 
displacernent or free surface displacement set and the complement of the total 
surface displacement set, respectively. The structural dynamic equation set 
with applied fluid pressure loading is, in partitioned form 
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(20) 



where the structural mass and stiffness matrices are expressed in parti- 
tioned form in accordance with the surface reference displacement set and 
the complement of the total surface displacement set, respectively. If the 
reference subsets of pressure, TP^, and displacement, XJ£, should correspond 
to a free fluid surface rather than a structural interface, then the structural 
mass partitions Mff, Mgf, Mfg would be null. In addition, the stiffness par- 
titions Kf£, Kfg, Kgf would correspond to contributions due to the surface 
gravitational potential and possibly ullage pressure fluctuation. In many 
cases the free surface and ullage stiffness are negligible relative to the struc- 
tural stiffness. K g, and may therefore be neglected. The development of 
the incompressible hydroelastic equations presented below treats the most 
general case implied in equation (20), with the above special cases illustrated 
at the completion of the derivationo 

From equation (19), the internal pressures are related to the surface 
pressures by 


P. 
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L"^[L..L. ] 
11 '• if IS-" 



( 21 ) 


and the reduced fluid dynamic equation set, in terms of surface quantities 
only, is 
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inertance matrix is singular since an incompressible fluid under 
umform pressure does not deform. This singularity is assured in the indi- 
vidual finite element inertance matrices by the use of admissible pressure 
distributions; it is a necessary condition for compatibility. 

For a fluid represented by discrete surface pressures the normalized 
umform pressure state is 


( 23 ) 


Under such loading the surface normal accelerations must be n ull and the 
necessary property of the reduced inertance matrix is 



( 24 ) 


■^e uniform pressure singularity pertains as well to the full pressure set with 
the members of Pj also unity; the reduction to the surface pressure set is 
made in this derivation mainly to illustrate the first basic property, of geo- 
rnetric dependence of internal pressure. For axisymmetric fluid elements 
with generalized pressure variables corresponding to amplitudes of circum- 
ferential harmonics, this discussion pertains to the zero«i harmonic, only, 
which contains the uniform pressure singularity. The singularity is removable 
by introduction of the concept of pressure deviation in which the deviations 
from the reference uniform pressure state of value, Pf, are 



( 25 ) 


Thus the relationship between pressures and pressure deviations may be ex- 
pressed by the transformation 
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( 26 ) 


when represents a single discrete reference pressure. If P^ represents a 
pressure subset comprised of aU free surface pressures (to be nulled in the 
case of negligible surface gravitational potential and ullage pressure strain 
energy), the constraint relationship is somewhat different; consideration of this 
latter case, however, is reserved for a later part of the present discussion. 
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Application of the above transformation on the reduced fluid equation 
set in equation (22), noting equation (24), yields 


II rr. .. T :: 

h, “s 


(27a) 


t I T* •• 'T' .. 

L P = - a/ Ua - A U 
ss s fs f ss s 


(27b) 


Suppose now that the displacement is directed normal to the reference sur- 
face; the area coupling partitions A^ and aT^ are therefore null, simplifying 
the above expressions. It should be noted that this restriction is ma.de for 
clarity of the present development and the more general case is derivable wicn 
some additional algebraic complication. Solution for the pressure deviations 
in equation (27b) and substitution of the result in equation (27a) yields the 
pressure deviation recovery relationship and the displacement recovery re- 
lationship, respectively. 


I ' 1 T “ 

P = - L U 

s ss ss s 


(28a) 


ri U 


(28b) 


with 


-T ' ' -1 T 

A,/ L, L A^^^ 
ff fs ss ss 


(28c) 


The latter result is equivalent to imposing a kinematic constraint on outward 
normal surface flow (incompressibility) 


a"^ U 

s s s 


I gi|a;^ uj 


(29a) 
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J I ss s 
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(29b) 
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the companion pressure (compatibility) constraint expression is 
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( 30 ) 


Physically the matrix G must consist of a row matrix with unit entries when 
Pf represents a single reference pressure. It is noted now that the use of an 
identity matrix in equation (28) and equation (30), implying that Pf and Uf may 
comprise pressure and acceleration subsets rather than individual quantities, 
is deliberate; the case in which a free surface pressure subset is null is 
therefore covered by the subsequent development. 

Substitution of equations (29) and (30) into equation (22a) yields 
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The lower partition merely represents the result already obtained in equation 
y u upper partition may be interpreted in two ways which consists of 

(1) the general case in which a single reference pressure, Pf, is chosen to 
express the uniform pressure state and (2) the special case in which a set of 
ree surface pressures, P£, must be zero. For the first general case, the 
singularity condition is expressed as 

^ff ■ ^fs Ss'^ ^sf = ° (31b) 


The constraints presented in equations (29) and (30) are now applied to the 
structural dynamic equation set equation (20), noting equations (28) and (31), 
resmting in the symmetric kinematic equation set 


(M" + M.) U + K'’‘ U = 0 
s IS s s 


where the fluid mass matrix is 
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and the constrained structural mass and stiffness matrices are 

m" = + m .r + m 

s ff fs sf ss 


and 


k'" = r^K„r 

S it 


r^K 


■fs 


+ K X 
sf 


+ K 


ss 


(32c) 


(32d) 


Complete displacement recovery is obtained through equation (28b) and pres- 
sure deviation recovery is obtained through equation (28a)o Surface pressure 
recovery is achieved by combining the upper partition equation set in equation 
(20) with equation (28b); thus the surface pressure recovery equation set con- 
sists of 


Pj = a;/ (MjjP + 0^ + a;/ (Kj^r + 


and equations (28a) and (30). 

The symmetric kinematic formulation developed above is useful in 
hydroelastic analyses for which either a fluid free surface is not present or 
free surface gravitational potential and/or ullage stiffness is significant. In 
most practical analyses involving tanks partially filled with fluid, the free 
surface strain energy is insignificant relative to the structural energy. In 
such cases low frequency slosh dynamics is approximated with rigid structure 
and flexible structure /fluid interaction dynamics is approximated with zero 
free surface pressure. In the latter case, the reference pressure set, Pp 
consists of all free surface pressures (set to zero) and the reference surface 
displacement set Uf consists of all free surface displacements. It is noted 
again as in the opening of this section that when a free surface is present Mff 
and Mfg are null; and when surface strain energy is insignificant, Kff and 
are null. 


DISPLACEMENT SET REDUCTION FOR THE 
HYDROELASTIC PROBLEM 


General Considerations 


For typical launch vehicle propellant tank models, the structural grid 
displacement set, Ug, may be fairly large (in excess of 1000 degrees of 
freedom) and the fluid mass matrix, equation (32b), is typically full. Reduc- 
tion to a much smaller dynamic set of variables is therefore very desirable 
for computational economy in free vibration analysis. The fullness of the 
mass matrix is peculiar to the hydroelastic problem and therefore special 
care must be taken in the reduction process. 
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Consider the reduction transformation 


U 

s 




(34) 


formed by static condensation of the stiffness matrix (Guyan reduction) or by 
some other process. Application of this transformation to the hydroelastic 
equation set, equation (32), with the constituent matrices explicitly defined 
yields the expressions for reduced mass and stiffness matrices . The most 
critical operations with respect to computational economy are those used in 
formation of the reduced fluid mass matrix 


(M^) 

reduced 




-1 T 

L ^ (A-^ T J 
ss ss sd 


(35) 


If the fluid mass matrix in the structural grid set, equation (32b), is expli- 
citly formed, a symmetric reduction process with T j may be quite time 
constiming and expensive due to the fullness of the omginal mass matrix. 
However, by first forming the reduced generalized area matrix aTt , as 
implied above in equation (35), the large, full fluid mass matrix neecf 
not be explicitly formed; the reduced mass matrix is calculated only, re- 
sulting in substantial computational economy. 


Harmonic Reduction of Geometrically 
Axisymmetric Structures 


The general category of structures of interest in propellant tank 
hydroelastic analysis consists of configurations which possess, for the most 
part, geometrically axisymmetric fluid cavities. The tank structure may 
have a variety of structural asymmetries such as circumferential thickness 
variation in both shell and rings, discrete longitudinal stiffeners and asym- 
metric supports. The finite element description of such a structure typi- 
cally must have a fine and fairly uniform nodal mesh distribution for adequate 
description of the structural behavior. The use of Guyan reduction in such 
cases (when the grid-set displacement degrees of freedom are on the order 
of thousands) may be inefficient and result in an inaccurate dynamic des- 
cription when extreme coordinate reduction is used. An alternative scheme, 
which expresses the circumferential displacement distribution in terms of a 
chosen set of harmonics (harmonic reduction), alleviates the difficulties 
encountered with Guyan reduction. This is especially true when only a few 
harmonic shapes appear sufficient to describe the anticipated dynamic be- 
havior of interest. 


Harmonic reduction of a discrete structural grid is accomplished by 
use of the geometric transformation (see ref. 9 for details) 


U 

s 




(36) 
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where Ug corresponds to the physical grid degrees of freedom to be trans- 
formed, and Ujj corresponds to the harmonic degrees of freedom (plus any 
discrete degrees of freedom not transformed). The transformation or 
constraint matrix is composed of the appropriate sinusoidal functions 
evaluated at the discrete variable locations. It should be noted that discrete 
displacements expressed in terms of cylindrical or spherical reference 
frames are most convenient for this procedure. For a typical shell struc- 
ture with JxK grid points such that there are J meridional rows and K cir- 
cumferential points in a row, the grid set has typically 6 x JxK degrees of 
freedom and the matrix semibandwidths are 6xK (assuming K < J). Appli- 
cation of the harmonic transformation as a reduction scheme, where the 
number of harmonics N is much less than K, results in a IL set of 6 x JxN 
generalized coordinates with matrix semibandwidths of 6xN. If N « K 
harmonic reduction represents a radical reduction in the number of degrees 
of freedom as well as matrix bandwidth. Further reduction of the system 
description is possible by a. small Guyan reduction made by choosing the 
generalized rotation degrees of freedom and tangential degrees of freedom 
as members of the omitted set of displacements. In such a case the analysis 
set consists of JxN degrees of freedom. This represents a radical reduction 
in degrees of freedom by a factor of (N/6K), without a costly large matrix 
decomposition typical of Guyan reduction. 


NUMERICAL STUDIES 


The new incompressible hydroelastic formulation and harmonic re- 
duction have been implemented in NASTRAN and verified and demonstrated 
on a number of problems. The problems fall into two categories, namely, 
analytical verification problems for which exact solutions are known, and 
demonstration problems for which experimental data are available. The 
1 /8-scale Space Shuttle external tank is included in the second category. The 
fluid idealizations utilized in the hydroelastic problems are based on existing 
elements of revolution in NASTRAN, 

The present hydroelastic analysis employs NASTRAN structural as 
well as fluid elements and provides for a description of dynamics of axi- 
symmetrically configured fluids in terms of circumferential harmonic pres- 
sure distributions. The distribution of pressure is typically 


P(r., e^.zj) = 


N 

k=l 


Ip (r., z.)cos k G. 
rk i’ i' 1 



z. )sin k G. 
1 1 


( 37 ) 
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The fluid containing structure is described in terms of discrete physical 
displacements so that the structural representation is not limited to struc- 
turally axisymmetric containers. Coupling of harmonic pressure distribu- 
tions with discrete structural displacements in this formulation is not strictly 
consistent; moreover, in many cases it is inefficient. When a structure 
described in terms of discrete displacements is coupled with a fluid described 
in terms of circumferential harmonics, inconsistencies may arise if too few 
pressure harmonics are utilized; structural deformation shapes associated 
with higher harmonics not included in the fluid representation will reflect a 
lack of fluid inertia loading. Alternatively, when the discrete structural 
grid is too coarse to accurately describe the highest harmonic pressure dis- 
tributions, large errors in the mode shapes associated with higher harmonics 
will be present. A consistent grid representation is realized by use of har- 
monic reduction when the number of structural harmonics coincides with the 
number of fluid pressure harmonics. This provides additional motivation 
for the use of harmonic reduction which is peculiar to use of the NASTRAN 
fluid elements. When harmonic reduction is not utilized, special care must 
be taken to utilize fluid harmonic and discrete structural descriptions of 
equivalent complexity. 


Analytical Verification Problems 


Spherical Cap . Harmonic reduction was first demonstrated on a 
spherical cap structure of uniform thickness to radius ratio, h/R = 0.05. 

The base of the cap, 60 degrees from the pole, is taken as rigidly fixed.* 

The material properties are: elastic modulus E = lO^, a Poisson's ratio 
V = 0. 3, and mass density = 0. 05. 

A nodal grid, consisting of 20 circumferential divisions in a semicircle 
and 10 meridional divisions, was chosen resulting in a structural model with 
1266 DOF (degrees of freedom). Three circumferential harmonics (0, 1, 2) 
were chosen for harmonic reduction. The apex node was left in terms of 
rectangular coordinates since the polar degrees of freedom have no meaning 
at this node. After application of the fixed-base boundary condition and 
symmetric kinematic constraints at the pole, and a small Guyan reduction, 
an analysis set of 72 outward-normal and meridional generalized harmonic 
displacements resulted. At this point, all natural frequencies and the first 
15 modes were calculated. Circumferential harmonics were uncoupled 
because of the axisymmetry in shell thickness. 

The results of the above approach were compared to results based on 
various Guyan reduction strategies and to "exact" results based on the 
STARS-II program (ref. 11). A comparison of computed natural frequencies 
( table I ) indicate that the overall accuracy of the 72 -DOF harmonic reduction 
representation is better than the 190-DOF Guyan reduction representation. 

The computation time associated with eigenvalue analysis of the harmonic 
analysis set is much less than that associated with the Guyan reduction 
analysis set. Central processing unit (CPU) times were 238 and 531 seconds, 
for the harmonic and Guyan reduction representations, respectively. This is 
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attributed to elimination of a large-scale matrix decomposition, character- 
istic of Guyan reduction, and to the fact that much fewer degrees of freedom 
are required for comparable accuracy (Sogo, harmonic 72-DOF for 5% accu- 
racy versus Guyan 190-DOF for 12% accuracy. 

Fluid-Filled Hemisphere . The first problem for hydroelastic verifica- 
tion consists of an open hemispherical container filled with fluid. The con- 
tainer is massless and follows the artificial structural law 


p = (38) 

where u^ is the local radial displacement; the exact free -vibration solution 
is known (ref, 12). The finite element model of the fluid and container is 
illustrated in figure 1, A diagonal structural stiffness matrix with entries 

K.. = aA. (39) 

11 1 


results from the artificial structural law in which Ai is the area associated 
with the ”ith" radial degree of freedom. The fluid model is expressed in 
terms of the circumferential pressure harmonics n = 0,2,4 and the struc- 
tural surface and free surface grids are reduced by harmonic reduction ac- 
cordingly. The fluid mass matrix is expressed in terms of a 21-DOF 
analysis set of structural radial displacements at seven meridional locations. 

A comparison of exact and finite element nondimensional natural fre- 
quencies is presented in table II and comparisons of selected modal dis- 
placement distributions are presented in figure 2. The finite element results 
are in very good agreement with the exact solution, with the level of accu- 
racy decreasing with modal complexity as expected. 

Fluid -Filled' Cylinder . Another hydroelastic verification problem 
consists of the fluid-filled, circular cylindrical shell illustrated in figure 3. 
The shell structure is taken as one with bending as well as membrane stiff- 
ness. The geometric properties of the shell consist of a cylinder with 
length/radius ratio (i/R) = 2 and a thickness /radius ratio (h/R) = 0.01, In 
addition, the ratio of fluid to structure density {Pr/P ) = 1/3 and the struc- 
tural material Poisson ratio (v) = 0,3. An exact nydroelastic modal solution 
is known for an infinitely long cylinder (refs. 13, 14) which holds for the 
present problem when the structure is subjected to the boundary conditions 


u = M = N = N =0 (shear diaphragm for 
^ ® ® ^ z = I r = R) 


= 0 (free surface) for 7 .- 1 , r £ R 


(40a) 

(40b) 
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9u^ 9M 

9^ ^ '^z ^ ~8i“ "" ^ ° (symmetry) for (40c) 

z = 0, r = R 


= 0 (fixed bottom) for z = 0, r :£ R 


The finite element models of the quarter shell (-0<z<i, 0<6< 

90 deg) and fluid were constructed taking advantage of symmetry. The 
structural grid for the quarter shell consists of 726 DOF (11 meridional 
nodal rows, 11 circumferential nodal columns) and the fluid grid consists of 
165 DOF (55 nodes of rotation, circumferential harmonics n = 0,2,4). A 
30-DOF set was obtained using (1) a harmonic transformation retaining har- 
monics n = 0,2,4, (2) the application of single-point constraints to enforce 
boundary conditions, and (3) a small Guyan reduction retaining only radial 
displacements. 

All 30 natural frequencies and 25 mode shapes with and without the 
fluid included were calculated. Frequency spectra for the empty and fluid- 
filled shells are presented in figure 4, illustrating generally excellent com- 
parison between finite element and exact results in both cases. 


A characteristic of the present formulation, which is as significant as 
numerical accuracy, is computational economy. On the IBM 370/165 com- 
puter the total solution time for the empty cylinder was about 2 CPU minutes 
for the fluid-filled cylinder an additional CPU minute was required to form 
the fluid mass and pressure recovery matrices. 


Comparisons with Experimental Data 


Liquid-Filled Cylinders Under Static Pressurization . A detailed ex- 
perimental study of the dynamics of structurally axisymmetric and asym- 
metric circular cylinders under various water fill and static pressurization 
conditions has been conducted at NASA Langley Research Center by Mr. 
Robert Herr. Data resulting from these tests (unpublished) are quite com- 
plete and provide an excellent basis for analysis /test correlation studies. 
The test articles are aluminum cylinders with mean radius of 2 5.4 cm 
(10 in. ) and height of 50. 8 cm (20 in. ). The cylinder walls are welded at 
the top and bottom to heavy aluminum plates. The axisymmetric test article 
has a cylinder wall thickness of 0. 081 cm (0. 032 in. ) and the asymmetric 
test article has a wall thickness variation around the circumference of 0, 051 
0. 102 cm (0. 020-0. 040 in. ) according to the equation 

= 0. 75 + 0. 25 cos 0 (41 ) 

max. 
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A structural model for the half-cylinder (0 deg - 0 £180 deg) taking 
advantage of the single axis of symmetry was assembled with sufficiently 
fine grid to simulate the circumferential harmonic shapes up to n = 15, which 
were known a priori from the experimental results to dominate in the lowest 
frequency modes. The structural grid set consists of 2046 DOF resulting 
from 11 axial grid rows and 31 circumferential grid columns all evenly dis- 
tributed, The fluid representation chosen at the half -filled condition is 
illustrated in figure 5. The 480-DOF (pressures) fluid model results from 
30 fluid-grid locations of revolution expressed in terms of the circum- 
ferential harmonics n = 0 to 15, with a sufficiently fine grid near the struc- 
tural wall to simulate the sharp pressure gradients occurring in the higher 
harmonics. Since modes with significant harmonic content up to n = 15 were 
of interest, harmonic reduction was not utilized here; computational economy 
would not be improved by the harmonic transformation and thus only a Guyan 
reduction was utilize do 

A series of cases including symmetric and asymmetric cylinders in 
the half-filled and empty configurations were studied. In addition, the ef- 
fects of static pressurization were included by use of differential stiffness 
capability in NASTRAN. 

The first cases studied pertained to the cylinder of uniform thickness. 
The empty cylinder was first considered with an assumed axial plane of 
symmetry at z = 25. 4 cm such that only m = 1, 3, 5 modes would be calcu- 
lated. The grid set of 1116 DOF, consisting of nodes below z = 25.4 cm, 
was reduced by Guyan reduction to an analysis set of 276 radial DOF with the 
lower end completely fixed (clamped). The unpressurized and pressurized 
gji©!!, frequency spectra of m — 1, n — 4 modes are illustrated in figure 6 
along with the test results. The calculated frequency spectrum was higher 
than the experimental frequency spectrum for both unpressurized and pres- 
surized conditions. A series of modifications of the structural model to 
reconcile the differences in results were considered. It was finally con- 
cluded that axial flexibility idealized as an axially free condition in the cylin- 
der/plate weld, provided the proper correction. Incorporation of the relaxed 
boundary conditions 


n = = N =0 for z = 0, 50. 8 cm (42) 

r dz z 

resulted in extremely accurate frequency spectra for the empty cylinder as 
illustrated in figure 6. 

The half-filled condition was then considered. The cylinder structure 
in this case does not have a dynamic plane of symmetry at z = 25.4 cm; the 
lower portion (z < 25.4 cm) is loaded by fluid structural inertia whereas the 
upper portion (z > 25.4 cm) is loaded only by the structural inertia. This 
provides motivation for Guyan reduction with all degrees of freedom at and 
above z =25.4 cm omitted (not including the supported degrees of freedom). 
A Guyan reduction on the structure and fluid was then performed resulting in 
an analysis set consisting of 248 radial DOF consistent with the fluid mass 
matrix. Hydroelastic modes, based on the clamped and modified-clamped 
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end conditions [eq. (40)], were then calculatedo The m = 1, n >4 modal 
frequency spectra, illustrated in figure 6, show that the representation with 
the relaxed end conditions is quite accurate, as concluded in the empty case* 

The unsymmetric cylinder structural model consists of the same grid 
set as in the case of the axisymmetric cylinder, and the Guyan reduction 
discussed above was utilized* The hydroelastic study of this cylinder was 
lirnited to the half-filled condition with the*'realistic”edge condition applied* 
The m = 1 mode shapes illustrated in figure 7 are in very good agreement 
with the unpressurized and pressurized test results, as are the modal fre- 
quencies. 

Computation times for the cylinder study were moderate since harmonic 
reduction, which was not appropriate, was not utilized. In all cases con- 
sidered, all 248 eigenvalues and 25 eigenvectors were calculated* Compu- 
tation time for the empty axisymmetric cylinder was 509 CPU sec. Prepa- 
ration of fluid matrix data req\oired 97 CPU sec and computation of hydro - 
elastic modes required 1,193 CPU sec. The increased CPU time required 
in this case is predominantly due to the increased structural grid set size 
of the fluid -filled cases; the increase in Guyan reduction time for systems 
of equivalent matrix bandwidth is proportional to the increase in grid set 
degrees of freedom. Computation times for the unsymmetric cylinder were 
similar to those required for the axisymmetric cylinder* 

The 1 /8-Scale Space Shuttle External Tank* An investigation of the dy- 
namic s^fanr78^^^scaIeSpane~ShuUle"’ex^^ in a free -free supported 

condition is in progress at NASA Langley Research Center. The 1 /8-scale 
external tank consists of two separate propellant tanks connected by a cylin- 
drical section. Although the fluid/structure interface is axisymmetric, the 
tank structure contains thickness and stiffener asymmetries* The finite 
element hydroelastic model for half the structure, taking advantage of the 
single axis of symmetry, is described in detail in references 9 and 15* It 
consists of a grid set of 348 pressure DOF and 2058 structural DOF (and 
768 harmonic structural DOF to be used in harmonic reduction). The 
structural grid deformed in the fundamental bending mode is illustrated in 
figure 8. Harmonics n = 0, 1,2,3 were chosen to describe asymmetric dy- 
namics with the pitch plane taken as an axis of symmetry. The analysis set 
of displacements resulting from a combination of harmonic and Guyan reduc- 
tions consists of 128 harmonic DOF associated with outward normal motion 
of the tank wall. 

Three conditions have been studied consisting of nearly full, interme- 
diate and empty propellant fill conditions. For each of the fill conditions, 

128 natural frequencies and 25 mode shapes and modal pressure distribu- 
tions were calculated with very good computational efficiency. About 20 
CPU minutes per liquid level on the IBM 370/165 computer was required to 
perform the entire analysis including matrix assembly, reduction and modal 
analysis. In previous attempts to study the dynamics of the same finite 
element representation with the standard unsymme'tric NAS TRAN hydro elas- 
tic analysis, computation times were about 52 CPU minutes for only one 
natural frequency and mode shape from a 412-DOF analysis set (ref. 15), 
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Excellent agreement between analysis and experimental frequencies 
occurred in the first axial mode but poor agreement occurred in the bending 
modes. A thorough investigation of the fluid model revealed good consis- 
tency in the idealization. As a result, the source of the discrepancy is be- 
lieved to be in the finite element representation of the structure which was 
constructed prior to the present work. 


CONCLUDING REMARKS 


Symmetric finite element matrix formulations for compressible and 
incompressible hydroelasticity have been developed on the basis of Toupin's 
complementary formulation of classical mechanics. The incompressible 
formulation applicable in propellant tank hydroelastic analysis has been 
implemented in NASTRAN to replace an inefficient \ansymmetric matrix for- 
mulation. The new technique which utilizes existing fluid and structural 
finite elements has been verified and demonstrated to be accurate and effi- 
cient. 


The fluid representation in the incompressible case reduces to a sym- 
metric fluid mass matrix described in terms of surface deformation only 
upon recognition of a singularity in the fluid inertance matrix. The singu- 
larity describes a physically necessary compatibility condition in that it 
assures that the incompressible idealization will not move under uniform 
pressure. Moreover, the singularity defines a kinematic constraint which is 
applied to the structural idealization when the fluid is completely bounded by 
a structural interface and when free surface ullage and/or gravitational 
stiffness are significant. The fluid mass matrix is added directly to the 
structural mass matrix, forming a symmetric set of hydroelastic equations 
in terms of structural displacements. Modal hydroelastic analysis is per- 
formed with the same efficiency as in the case of a non-fluid-filled structure, 
since no additional degrees of freedom are required for the fluid (other than 
free surface displacements when necessary). 

The efficiency of the new hydroelastic analysis technique has been en- 
hanced for both fluid and structure by introduction of harmonic reduction, 
applicable to geometrically axisymmetric structures, as an alternative to 
Guyan reduction. When the number of harmonics utilized is much less than 
the number of discrete nodes about a circumference, overall matrix size 
and bandwidth are significantly reduced. 

The formulation has been verified by comparison with exact analytical 
results for a fluid-filled hemispherical container and a fluid-filled circular 
cylindrical shell. In all cases, excellent correlation was exhibited as well 
as very good computational efficiency. In addition, the analysis /test cor- 
relation study on symmetric and unsymmetric circular cylindrical shells 
under various fluid-fill conditions is considered very good. 

Analysis /test discrepancies on the l/8-scale external tank model for 
the space shuttle have not yet been resolved. The efficiency of the current 
analysis, however, is very encouraging based upon comparison of 
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computation times between the present and the standard unsymmetric 
NASTRAN hydroelastic formulations. 

The more general kinetic formulation which includes fluid compressi- 
bility has yet to be investigated in detail and applied. Typical applications 
include underwater explosion and acoustic analysis in general. In addition, 
a set of polyhedral complementary fluid finite elements should be incorporated 
to allow for modeling of general fluid configurations. 
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TABLE I. - 60-DEG SPHERICAL CAP - COMPARISON OF 
MODAL FREQUENCIES 



•in ^ 

Error = 1 

exact 


TABLE II. - FLUID IN A HEMISPHERICAL CONTAINER 
NATURAL FREQUENCY COMPARISONS 



1 










FLUID GRID OF REVOLUTION STRUCTURAL GRID FOR 

QUARTER HEMISPHERE 



FIGURE 1. - FLUID AND HEMISPHERICAL 
CONTAINER MODEL 





• DENOTES FINITE ELEMENT MODE SHAPE 


FIGURE 2. - HEMISPHERE HYDROELASTIC 
MODE SHAPES 
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FIGURE 3. - FLUID-FILLED CIRCULAR CYLINDRICAL SHELL 



AXIAL WAVELENGTH PARAMETER, 8/mR AXIAL WAVELENGTH PARAMETER, 8/mR 

■“ "EXACT ' THEORY 
• FINITE ELEMENT SOLN 


FIGURE 4. - CIRCULAR CYLINDRICAL SHELL FREQUENCY SPECTRA 





p = 0 ON FREE SURFACE 


z 



I 1 — I — ^ r 

•' = 0 r = 25.4 cm 


FIGURE 5. - HALF-FILLED CYLINDER FLUID IDEALIZATION 

A TEST 

□ ANALYSIS. FULLY CLAMPED ENDS 
O ANALYSIS, RELAXED END RESTRAINTS 


200 ' ' ' 1 1 

46 8 10 12 

n 



FIGURE 6. - AXISYMMETRIC CYLINDER FREQUENCY 
SPECTRA (m = 1 MODES) 
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0. 102 cm 


0.051 ( 


THICKNESS DISTRIBUTION EQ- 4 1 
0 ° 180 ' 



z = 20. 32 cm 




(a) NO PRESSURIZATION, p=0 




MODE 2: ~ ^TEST”^^^ 


lb) WITH PRESSURIZATION, = 5. 5 x 10^ N/M* 


FIGURE 7. - UNSYMMETRIC CYLINDER HYDROELASTIC MODES 



FIGURE 8. - 1 /8-SCALE SPACE SHUTTLE EXTERNAL TANK 
FUNDAMENTAL HYDROELASTIC BENDING MODE 
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